Mushroom Data Analysis

Índice

1. Introducción
   
2. Cargar librerías
   
3. Cargar datos
   
4. Visualización y preprocesamiento
   4.1 Imputación de valores nulos
   4.2 Codificación de variables categóricas
   4.3 Escalado de variables numéricas
5. Análisis de datos
   5.1 Análisis aprendizaje supervisado
   5.2 Análisis aprendizaje no supervisado
   5.2.1 K-means
   5.2.2 Clustering jerárquico
   

Introduction

Los datos provienen de la página web Kaggle. El dataset contiene información sobre setas comestibles y venenosas. El dataset contiene 8124 observaciones y 24 variables. Las variables con las que trabajaremos son las siguientes:

1. class: edible=e, poisonous=p
2. cap-diameter: continuous
3. cap-shape: bell=b, conical=c, convex=x, flat=f, knobbed=k, sunken=s
4. cap-surface: fibrous=f, grooves=g, scaly=y, smooth=s
5. cap-color: brown=n, buff=b, cinnamon=c, gray=g, green=r, pink=p, purple=u, red=e, white=w, yellow=y
6. bruises: bruises=t, no=f
7. odor: almond=a, anise=l, creosote=c, fishy=y, foul=f, musty=m, none=n, pungent=p, spicy=s
8. gill-attachment: attached=a, descending=d, free=f, notched=n
9. gill-spacing: close=c, crowded=w, distant=d
10. gill-size: broad=b, narrow=n
11. gill-color: black=k, brown=n, buff=b, chocolate=h, gray=g, green=r, orange=o, pink=p, purple=u, red=e, white=w, yellow=y
12. stalk-shape: enlarging=e, tapering=t
13. stalk-root: bulbous=b, club=c, cup=u, equal=e, rhizomorphs=z, rooted=r, missing=?
14. stalk-surface-above-ring: fibrous=f, scaly=y, silky=k, smooth=s
15. stalk-surface-below-ring: fibrous=f, scaly=y, silky=k, smooth=s
16. stalk-color-above-ring: brown=n, buff=b, cinnamon=c, gray=g, orange=o, pink=p, red=e, white=w, yellow=y
17. stalk-color-below-ring: brown=n, buff=b, cinnamon=c, gray=g, orange=o, pink=p, red=e, white=w, yellow=y
18. veil-type: partial=p, universal=u
19. veil-color: brown=n, orange=o, white=w, yellow=y
20. ring-number: none=n, one=o, two=t
21. ring-type: cobwebby=c, evanescent=e, flaring=f, large=l, none=n, pendant=p, sheathing=s, zone=z
22. spore-print-color: black=k, brown=n, buff=b, chocolate=h, green=r, orange=o, purple=u, white=w, yellow=y
23. population: abundant=a, clustered=c, numerous=n, scattered=s, several=v, solitary=y
24. habitat: grasses=g, leaves=l, meadows=m, paths=p, urban=u, waste=w, woods=d

Cargar librerías

install.packages("caret")

trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/caret_6.0-93.tgz'
Content type 'application/x-gzip' length 3578325 bytes (3.4 MB)
==================================================
downloaded 3.4 MB

The downloaded binary packages are in
    /var/folders/v0/47swnvh93cl0_bw1dw8tl_0h0000gn/T//RtmpBAcuFX/downloaded_packages

install.packages("tidyverse")

trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/tidyverse_1.3.2.tgz'
Content type 'application/x-gzip' length 425892 bytes (415 KB)
==================================================
downloaded 415 KB

The downloaded binary packages are in
    /var/folders/v0/47swnvh93cl0_bw1dw8tl_0h0000gn/T//RtmpBAcuFX/downloaded_packages

install.packages("plotly")

trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/plotly_4.10.1.tgz'
Content type 'application/x-gzip' length 3179590 bytes (3.0 MB)
==================================================
downloaded 3.0 MB

The downloaded binary packages are in
    /var/folders/v0/47swnvh93cl0_bw1dw8tl_0h0000gn/T//RtmpBAcuFX/downloaded_packages

install.packages("dplyr")

trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/dplyr_1.0.10.tgz'
Content type 'application/x-gzip' length 1327807 bytes (1.3 MB)
==================================================
downloaded 1.3 MB

The downloaded binary packages are in
    /var/folders/v0/47swnvh93cl0_bw1dw8tl_0h0000gn/T//RtmpBAcuFX/downloaded_packages

install.packages("factoextra")

trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/factoextra_1.0.7.tgz'
Content type 'application/x-gzip' length 415155 bytes (405 KB)
==================================================
downloaded 405 KB

The downloaded binary packages are in
    /var/folders/v0/47swnvh93cl0_bw1dw8tl_0h0000gn/T//RtmpBAcuFX/downloaded_packages

install.packages("dendextend")

trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/dendextend_1.16.0.tgz'
Content type 'application/x-gzip' length 3895351 bytes (3.7 MB)
==================================================
downloaded 3.7 MB

The downloaded binary packages are in
    /var/folders/v0/47swnvh93cl0_bw1dw8tl_0h0000gn/T//RtmpBAcuFX/downloaded_packages

library(caret)

Loading required package: ggplot2
Loading required package: lattice
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.1 
✔ readr   2.1.3      ✔ forcats 0.5.2 
✔ purrr   0.3.5      ── Conflicts ─────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
✖ purrr::lift()   masks caret::lift()

library(plotly)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(dplyr)
library(cluster)
library(factoextra)

Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

library(dendextend)

---------------------
Welcome to dendextend version 1.16.0
Type citation('dendextend') for how to cite the package.

Type browseVignettes(package = 'dendextend') for the package vignette.
The github page is: https://github.com/talgalili/dendextend/

Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
You may ask questions at stackoverflow, use the r and dendextend tags: 
     https://stackoverflow.com/questions/tagged/dendextend

    To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
---------------------


Attaching package: ‘dendextend’

The following object is masked from ‘package:stats’:

    cutree

Cargar datos

En primer lugar, realizamos la lectura del fichero csv, separándolo por “;” y mostramos las primeras 6 filas del fichero.

mushroom <- read.csv("./data/data.csv", sep = ";")
head(mushroom)

Echamos un vistazo a las características de los atributos del dataset. En el caso de las variables numéricas, se puede observar valores como el mínimo, máximo, media, desviación estándar, etc.

summary(mushroom)

    class            cap.diameter     cap.shape         cap.surface       
 Length:61069       Min.   : 0.380   Length:61069       Length:61069      
 Class :character   1st Qu.: 3.480   Class :character   Class :character  
 Mode  :character   Median : 5.860   Mode  :character   Mode  :character  
                    Mean   : 6.734                                        
                    3rd Qu.: 8.540                                        
                    Max.   :62.340                                        
  cap.color         does.bruise.or.bleed gill.attachment    gill.spacing      
 Length:61069       Length:61069         Length:61069       Length:61069      
 Class :character   Class :character     Class :character   Class :character  
 Mode  :character   Mode  :character     Mode  :character   Mode  :character  
                                                                              
                                                                              
                                                                              
  gill.color         stem.height       stem.width      stem.root         stem.surface      
 Length:61069       Min.   : 0.000   Min.   :  0.00   Length:61069       Length:61069      
 Class :character   1st Qu.: 4.640   1st Qu.:  5.21   Class :character   Class :character  
 Mode  :character   Median : 5.950   Median : 10.19   Mode  :character   Mode  :character  
                    Mean   : 6.582   Mean   : 12.15                                        
                    3rd Qu.: 7.740   3rd Qu.: 16.57                                        
                    Max.   :33.920   Max.   :103.91                                        
  stem.color         veil.type          veil.color          has.ring        
 Length:61069       Length:61069       Length:61069       Length:61069      
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
  ring.type         spore.print.color    habitat             season         
 Length:61069       Length:61069       Length:61069       Length:61069      
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            

Visualización y preprocesamiento

Debido a que el dataset contiene valores vacíos, los sustituiremos por valores NA para poder trabajar con ellos.

mushroom[mushroom == ""] <- NA

Comprobamos que se han sustituido correctamente los valores vacíos por NA.

str(mushroom)

'data.frame':   61069 obs. of  21 variables:
 $ class               : chr  "p" "p" "p" "p" ...
 $ cap.diameter        : num  15.3 16.6 14.1 14.2 14.6 ...
 $ cap.shape           : chr  "x" "x" "x" "f" ...
 $ cap.surface         : chr  "g" "g" "g" "h" ...
 $ cap.color           : chr  "o" "o" "o" "e" ...
 $ does.bruise.or.bleed: chr  "f" "f" "f" "f" ...
 $ gill.attachment     : chr  "e" "e" "e" "e" ...
 $ gill.spacing        : chr  NA NA NA NA ...
 $ gill.color          : chr  "w" "w" "w" "w" ...
 $ stem.height         : num  16.9 18 17.8 15.8 16.5 ...
 $ stem.width          : num  17.1 18.2 17.7 16 17.2 ...
 $ stem.root           : chr  "s" "s" "s" "s" ...
 $ stem.surface        : chr  "y" "y" "y" "y" ...
 $ stem.color          : chr  "w" "w" "w" "w" ...
 $ veil.type           : chr  "u" "u" "u" "u" ...
 $ veil.color          : chr  "w" "w" "w" "w" ...
 $ has.ring            : chr  "t" "t" "t" "t" ...
 $ ring.type           : chr  "g" "g" "g" "p" ...
 $ spore.print.color   : chr  NA NA NA NA ...
 $ habitat             : chr  "d" "d" "d" "d" ...
 $ season              : chr  "w" "u" "w" "w" ...

Comprobaremos la cantidad de valores nulos que hay en cada variable. Para ello, utilizaremos la función colSums(is.na(mushroom)), que nos devolverá la suma de valores nulos de cada variable. Esta información nos ayudará a decidir si será conveniente eliminar dichas variables o no. Podemos observar que existen 5 variables donde más del 50% de los valores son nulos. Estas son: stem.surface, veil.color, spore.print.color, stem.root, veil.type. Esta información nos ayudará a decidir si será conveniente eliminar dichas variables o no. En este caso, dichas variables serán eliminadas ya que, en caso de imputar los valores nulos, se perdería mucha información.

colSums(is.na(mushroom))

               class         cap.diameter            cap.shape          cap.surface 
                   0                    0                    0                14120 
           cap.color does.bruise.or.bleed      gill.attachment         gill.spacing 
                   0                    0                 9884                25063 
          gill.color          stem.height           stem.width            stem.root 
                   0                    0                    0                51538 
        stem.surface           stem.color            veil.type           veil.color 
               38124                    0                57892                53656 
            has.ring            ring.type    spore.print.color              habitat 
                   0                 2471                54715                    0 
              season 
                   0 

En este caso, decidimos eliminar aquellas variables con más del 50% de valores nulos ya que, en caso de imputar los valores nulos, se perdería mucha información. Para ello, en primer lugar, obtendremos el nombre de las columnas que queremos eliminar.

nacols <- colnames(mushroom)[colSums(is.na(mushroom)) > nrow(mushroom) / 2]
print(nacols)

[1] "stem.root"         "stem.surface"      "veil.type"         "veil.color"       
[5] "spore.print.color"

Eliminamos las variables no deseadas.

mushroom <- mushroom[, !names(mushroom) %in% nacols]

Comprobamos que se han eliminado las variables con más del 50% de valores nulos.

print(colnames(mushroom))

 [1] "class"                "cap.diameter"         "cap.shape"           
 [4] "cap.surface"          "cap.color"            "does.bruise.or.bleed"
 [7] "gill.attachment"      "gill.spacing"         "gill.color"          
[10] "stem.height"          "stem.width"           "stem.color"          
[13] "has.ring"             "ring.type"            "habitat"             
[16] "season"              

A continuación, separamos las variables numéricas de las categóricas.

colsnames <- colnames(mushroom)
numerical_features <- c("cap.diameter", "stem.height", "stem.width")
categorical_features <- colsnames[!colsnames %in% numerical_features]
print(categorical_features)

 [1] "class"                "cap.shape"            "cap.surface"         
 [4] "cap.color"            "does.bruise.or.bleed" "gill.attachment"     
 [7] "gill.spacing"         "gill.color"           "stem.color"          
[10] "has.ring"             "ring.type"            "habitat"             
[13] "season"              

print(numerical_features)

[1] "cap.diameter" "stem.height"  "stem.width"  

Visualizamos la distribución de las variables categóricas a través de histogramas. Observaremos los posibles valores de cada variable categórica junto con su frecuencia de aparición en el dataset.

for (i in categorical_features) {
  print(ggplot(mushroom, aes_string(x = i)) +
    geom_bar(fill = "#7fd6d9") +
    geom_text(stat = "count", aes(label = scales::percent(..count.. / nrow(mushroom)), vjust = -0.25)) +
    labs(x = i, y = "Percentage") +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)))
}

Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
Please use tidy evaluation ideoms with `aes()`

Tras visualizar los histogramas de cada variable, llegamos a las siguientes conclusiones:

• La variable dependiente “class” está balanceada, es decir, la frecuencia de aparición de “e” (edible) y “p” (poisonous) es similar, en una proporción de 44,5% y 55,5% respectivamente. Por lo tanto, no es necesario realizar un balanceo de la variable dependiente, permitiendo aplicar la medida de evaluación “accuracy” para evaluar el modelo.
• Los valores NA de los atributos restantes serán imputados con la moda de cada variable.

A continuación, visualizaremos la distribución de las variables numéricas a través de una gráfica generada con “featurePlot”. Esta requiere que la variable objetivo sea de tipo factor, por lo que hacemos la conversión.

mushroom$class <- as.factor(mushroom$class)
featurePlot(x = mushroom[, numerical_features], y = mushroom$class, plot = "strip")

Con la gráfica anterior, se puede observar que para valores altos de cap.diameter, stem.height y stem.width, la probabilidad de que la clase sea “e” es mayor que la de “p”.

Para decidir si es conveniente eliminar alguna variable, se puede utilizar la función “nearZeroVar”. Esta función devuelve un vector con los índices de las variables que tienen una varianza cercana a 0.

near_zero_col <- nearZeroVar(mushroom, saveMetrics = FALSE)
colnames(mushroom)[c(near_zero_col)]

[1] "ring.type"

Como conclusión, la variable “ring.type” podría ser eliminada ya que tiene una varianza cercana a 0. A continuación, visualizaremos la correlación existente con la variable dependiente “class”.

print(ggplot(mushroom, aes_string(x = "ring.type")) +
  geom_bar(aes(fill = class)))

Tras visualizar la gráfica anterior, se puede observar que la distribución de la variable dependiente “class” es similar en casi todos los valores de “ring.type”. Sin embargo, en el caso de “ring.type” = “z” y “ring.type” = “m”, la distribución de la variable dependiente “class” es diferente. Por lo tanto, se decide no eliminar la variable “ring.type”.

Imputación de valores nulos

Como se ha comentado anteriormente, los valores nulos de las variables categóricas se imputarán a través de la moda de cada variable. Esto es debido a que la función preProcess() de la librería caret no permite imputar valores nulos de variables categóricas.

for (i in categorical_features) {
  mushroom[, i][is.na(mushroom[, i])] <- names(which.max(table(mushroom[, i])))
}

Podemos comprobar que ya no existen valores nulos en las variables categóricas.

colSums(is.na(mushroom))

               class         cap.diameter            cap.shape          cap.surface 
                   0                    0                    0                    0 
           cap.color does.bruise.or.bleed      gill.attachment         gill.spacing 
                   0                    0                    0                    0 
          gill.color          stem.height           stem.width           stem.color 
                   0                    0                    0                    0 
            has.ring            ring.type              habitat               season 
                   0                    0                    0                    0 

Escalado de variables numéricas

Para escalar las variables numéricas, se utilizará la función preProcess() de la librería caret. Esta función devuelve un objeto de tipo “preProcess” que contiene la información necesaria para escalar las variables numéricas. A continuación, se escalarán las variables numéricas y se eliminarán las variables originales.

range_numeric <- preProcess(mushroom[, numerical_features], method = c("range"))
mushroom[, numerical_features] <- predict(range_numeric, newdata = mushroom[, numerical_features])
str(mushroom)

'data.frame':   61069 obs. of  16 variables:
 $ class               : Factor w/ 2 levels "e","p": 2 2 2 2 2 2 2 2 2 2 ...
 $ cap.diameter        : num  0.24 0.262 0.221 0.223 0.23 ...
 $ cap.shape           : chr  "x" "x" "x" "f" ...
 $ cap.surface         : chr  "g" "g" "g" "h" ...
 $ cap.color           : chr  "o" "o" "o" "e" ...
 $ does.bruise.or.bleed: chr  "f" "f" "f" "f" ...
 $ gill.attachment     : chr  "e" "e" "e" "e" ...
 $ gill.spacing        : chr  "c" "c" "c" "c" ...
 $ gill.color          : chr  "w" "w" "w" "w" ...
 $ stem.height         : num  0.5 0.53 0.525 0.465 0.487 ...
 $ stem.width          : num  0.164 0.175 0.171 0.154 0.166 ...
 $ stem.color          : chr  "w" "w" "w" "w" ...
 $ has.ring            : chr  "t" "t" "t" "t" ...
 $ ring.type           : chr  "g" "g" "g" "p" ...
 $ habitat             : chr  "d" "d" "d" "d" ...
 $ season              : chr  "w" "u" "w" "w" ...

Análisis de datos

Análisis aprendizaje supervisado

Regresión Logística

Primero se dividirá el dataset en dos conjuntos: uno de entrenamiento y otro de test. El conjunto de entrenamiento se utilizará para entrenar los modelos y el conjunto de test se utilizará para evaluarlos.

library(caTools)
set.seed(18)

split <- sample.split(mushroom$class, SplitRatio = 0.8)
training_set <- subset(mushroom, split == TRUE)
test_set <- subset(mushroom, split == FALSE)

table(training_set$class)

    e     p 
21745 27110 

table(test_set$class)

   e    p 
5436 6778 

(Redactar) La regresión logística permite predecir el resultado de una variable categórica en función de las variables independientes o predictoras.

rl_classiffier <- glm(class ~ ., family = binomial, data = training_set)
summary(rl_classiffier)

Call:
glm(formula = class ~ ., family = binomial, data = training_set)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.98113  -0.77757   0.00045   0.75561   2.99676  

Coefficients: (2 not defined because of singularities)
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)           -16.146622 310.494270  -0.052 0.958526    
cap.diameter           -3.097299   0.273141 -11.340  < 2e-16 ***
cap.shapec             -1.587320   0.089853 -17.666  < 2e-16 ***
cap.shapef             -1.548574   0.056840 -27.244  < 2e-16 ***
cap.shapeo             -0.505473   0.102060  -4.953 7.32e-07 ***
cap.shapep             -1.355581   0.077525 -17.486  < 2e-16 ***
cap.shapes             -1.914206   0.067428 -28.389  < 2e-16 ***
cap.shapex             -1.587070   0.052730 -30.098  < 2e-16 ***
cap.surfacee            0.777609   0.081079   9.591  < 2e-16 ***
cap.surfaceg           -0.469069   0.067815  -6.917 4.62e-12 ***
cap.surfaceh           -0.869361   0.063758 -13.635  < 2e-16 ***
cap.surfacei            1.878605   0.116592  16.113  < 2e-16 ***
cap.surfacek            3.205155   0.106739  30.028  < 2e-16 ***
cap.surfacel           -1.338110   0.100947 -13.256  < 2e-16 ***
cap.surfaces           -0.913831   0.059362 -15.394  < 2e-16 ***
cap.surfacet           -0.037125   0.050461  -0.736 0.461905    
cap.surfacew           -0.483820   0.084640  -5.716 1.09e-08 ***
cap.surfacey           -0.579375   0.064113  -9.037  < 2e-16 ***
cap.colore              1.986990   0.109951  18.072  < 2e-16 ***
cap.colorg              0.582318   0.108490   5.367 7.98e-08 ***
cap.colork              1.553011   0.131690  11.793  < 2e-16 ***
cap.colorl             -0.129820   0.140094  -0.927 0.354103    
cap.colorn              0.374666   0.099619   3.761 0.000169 ***
cap.coloro              1.627975   0.112745  14.439  < 2e-16 ***
cap.colorp              1.051060   0.127472   8.245  < 2e-16 ***
cap.colorr              2.925006   0.136340  21.454  < 2e-16 ***
cap.coloru              1.538956   0.121348  12.682  < 2e-16 ***
cap.colorw              0.861936   0.104891   8.217  < 2e-16 ***
cap.colory              0.615058   0.102936   5.975 2.30e-09 ***
does.bruise.or.bleedt  -0.147991   0.038250  -3.869 0.000109 ***
gill.attachmentd        0.644857   0.046979  13.726  < 2e-16 ***
gill.attachmente       -0.922859   0.055651 -16.583  < 2e-16 ***
gill.attachmentf        0.816396   0.136750   5.970 2.37e-09 ***
gill.attachmentp       -2.450165   0.060717 -40.354  < 2e-16 ***
gill.attachments        0.142043   0.052091   2.727 0.006395 ** 
gill.attachmentx        0.089014   0.043392   2.051 0.040231 *  
gill.spacingd          -0.425196   0.036863 -11.534  < 2e-16 ***
gill.spacingf                 NA         NA      NA       NA    
gill.colore             2.304077   0.156835  14.691  < 2e-16 ***
gill.colorf                   NA         NA      NA       NA    
gill.colorg             0.833125   0.124899   6.670 2.55e-11 ***
gill.colork             0.854961   0.136405   6.268 3.66e-10 ***
gill.colorn             1.350254   0.119303  11.318  < 2e-16 ***
gill.coloro             1.131719   0.123598   9.156  < 2e-16 ***
gill.colorp             0.848246   0.122841   6.905 5.01e-12 ***
gill.colorr             0.896493   0.144108   6.221 4.94e-10 ***
gill.coloru             1.319431   0.145603   9.062  < 2e-16 ***
gill.colorw             0.819802   0.115149   7.120 1.08e-12 ***
gill.colory             1.818426   0.119373  15.233  < 2e-16 ***
stem.height             3.580359   0.196103  18.258  < 2e-16 ***
stem.width             -0.455339   0.225084  -2.023 0.043076 *  
stem.colore            17.897727 310.494233   0.058 0.954033    
stem.colorf            35.399323 334.067154   0.106 0.915610    
stem.colorg            15.906809 310.494228   0.051 0.959142    
stem.colork            19.545176 310.494262   0.063 0.949807    
stem.colorl            15.591779 310.494286   0.050 0.959950    
stem.colorn            17.491009 310.494224   0.056 0.955077    
stem.coloro            16.096825 310.494231   0.052 0.958654    
stem.colorp            19.131727 310.494251   0.062 0.950868    
stem.colorr            17.537012 310.494253   0.056 0.954959    
stem.coloru            17.419638 310.494223   0.056 0.955260    
stem.colorw            16.082601 310.494224   0.052 0.958691    
stem.colory            17.437089 310.494226   0.056 0.955215    
has.ringt              -0.006793   0.050561  -0.134 0.893124    
ring.typef             -0.911592   0.084553 -10.781  < 2e-16 ***
ring.typeg             -0.499679   0.101534  -4.921 8.60e-07 ***
ring.typel             -0.167323   0.100958  -1.657 0.097449 .  
ring.typem            -20.187930 230.159347  -0.088 0.930105    
ring.typep              0.894990   0.103503   8.647  < 2e-16 ***
ring.typer             -0.492646   0.105898  -4.652 3.29e-06 ***
ring.typez             16.741684  78.753199   0.213 0.831651    
habitatg                0.548824   0.041040  13.373  < 2e-16 ***
habitath                0.199731   0.065306   3.058 0.002225 ** 
habitatl               -0.501569   0.054900  -9.136  < 2e-16 ***
habitatm                0.103008   0.067655   1.523 0.127870    
habitatp               17.002090 230.326776   0.074 0.941156    
habitatu              -17.054326 378.536784  -0.045 0.964065    
habitatw              -17.436003 215.700665  -0.081 0.935574    
seasons                -1.618439   0.069225 -23.379  < 2e-16 ***
seasonu                 0.181688   0.025335   7.171 7.42e-13 ***
seasonw                -1.277528   0.049110 -26.014  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 67137  on 48854  degrees of freedom
Residual deviance: 45160  on 48776  degrees of freedom
AIC: 45318

Number of Fisher Scoring iterations: 16

Una vez hemos aplicado la función glm(), obtenemos los valores residuales del modelo y los coeficientes de ajuste para cada una de las variables independientes. Además, obtenemos el p-value correspondiente a cada una de ellas. Las variables con dos o tres asteriscos aportan bastante relevancia como predictores; sin embargo, las variables que no poseen ninguno o uno son de menor relevancia.

A continuación, procedemos a predecir las clases del conjunto de entrenamiento y de validación. Tomamos como umbral 0,5, de forma que si la probabilidad queda por encima de dicho umbral es que es comestible, pero si queda por debajo es que es venenoso. Por último crearemos la matriz de confusión para valorar los resultados de las predicciones.

pred_train <- predict(rl_classiffier, newdata = training_set, type = "response")

Warning: prediction from a rank-deficient fit may be misleading

pred_train <- ifelse(pred_train > 0.5, "p", "e")
pred_train <- factor(pred_train, levels = c("e", "p"), labels = c("e", "p"))

confusion_m <- table(training_set$class, pred_train)
confusion_m

   pred_train
        e     p
  e 16489  5256
  p  5573 21537

accuracy <- sum(diag(confusion_m)) / sum(confusion_m)
accuracy

[1] 0.7783441

Una vez hemos visto el resultado que muestra la matriz de confusión, se observa que la predicción tiene una precisión del 77,8%, lo que puede resultar algo baja.

Repetimos el mismo proceso para el conjunto de test.

pred_test <- predict(rl_classiffier, newdata = test_set, type = "response")

Warning: prediction from a rank-deficient fit may be misleading

pred_test <- ifelse(pred_test > 0.5, "p", "e")
pred_test <- factor(pred_test, levels = c("e", "p"), labels = c("e", "p"))

confusion_m <- table(test_set$class, pred_test)
confusion_m

   pred_test
       e    p
  e 4087 1349
  p 1420 5358

accuracy_rl <- sum(diag(confusion_m)) / sum(confusion_m)
accuracy_rl

[1] 0.7732929

Obtenemos una precisión con gran similitud a la del conjunto de training, un 77,3%. Para finalizar la regresión logística, vamos a graficar la curva ROC, la cual nos aporta la visualización de la relación entre los falsos y verdaderos positivos.

library(ROCR)
pred_rl_roc <- prediction(as.numeric(pred_test), as.numeric(test_set$class))
perf_rl_roc <- performance(pred_rl_roc, "tpr", "fpr")
perf_rl_auc <- performance(pred_rl_roc, "auc")

print(perf_rl_auc@y.values[[1]])

[1] 0.7711691

plot(perf_rl_roc)

Observando la curva ROC resultante podemos comentar que se mantiene por encima de la diagonal, lo que es buena señal, pero se aproxima a ella, pudiendo haber proporcionado resultados mejores.

k-NN

Antes de nada, para poder aplicar knn, debemos transformar las variables categóricas en numéricas.

mushroom_num <- dummyVars(" ~ .", data = mushroom, fullRank = TRUE) %>% predict(mushroom)
mushroom_num <- as.data.frame(mushroom_num)

Visualizamos las variables categóricas codificadas y las dimensiones del dataset.

dim_mushroom <- dim(mushroom_num)
print(dim_mushroom)

[1] 61069    81

str(mushroom_num)

'data.frame':   61069 obs. of  81 variables:
 $ class.p              : num  1 1 1 1 1 1 1 1 1 1 ...
 $ cap.diameter         : num  0.24 0.262 0.221 0.223 0.23 ...
 $ cap.shapec           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.shapef           : num  0 0 0 1 0 0 1 0 1 1 ...
 $ cap.shapeo           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.shapep           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.shapes           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.shapex           : num  1 1 1 0 1 1 0 1 0 0 ...
 $ cap.surfacee         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfaceg         : num  1 1 1 0 0 1 0 0 1 1 ...
 $ cap.surfaceh         : num  0 0 0 1 1 0 1 1 0 0 ...
 $ cap.surfacei         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfacek         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfacel         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfaces         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfacet         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfacew         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.surfacey         : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colore           : num  0 0 0 1 0 0 0 1 0 1 ...
 $ cap.colorg           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colork           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colorl           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colorn           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.coloro           : num  1 1 1 0 1 1 1 0 1 0 ...
 $ cap.colorp           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colorr           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.coloru           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colorw           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ cap.colory           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ does.bruise.or.bleedt: num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.attachmentd     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.attachmente     : num  1 1 1 1 1 1 1 1 1 1 ...
 $ gill.attachmentf     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.attachmentp     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.attachments     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.attachmentx     : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.spacingd        : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.spacingf        : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colore          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colorf          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colorg          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colork          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colorn          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.coloro          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colorp          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colorr          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.coloru          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ gill.colorw          : num  1 1 1 1 1 1 1 1 1 1 ...
 $ gill.colory          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.height          : num  0.5 0.53 0.525 0.465 0.487 ...
 $ stem.width           : num  0.164 0.175 0.171 0.154 0.166 ...
 $ stem.colore          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorf          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorg          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colork          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorl          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorn          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.coloro          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorp          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorr          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.coloru          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ stem.colorw          : num  1 1 1 1 1 1 1 1 1 1 ...
 $ stem.colory          : num  0 0 0 0 0 0 0 0 0 0 ...
 $ has.ringt            : num  1 1 1 1 1 1 1 1 1 1 ...
 $ ring.typef           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ring.typeg           : num  1 1 1 0 0 0 1 0 0 0 ...
 $ ring.typel           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ring.typem           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ring.typep           : num  0 0 0 1 1 1 0 1 1 1 ...
 $ ring.typer           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ ring.typez           : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitatg             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitath             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitatl             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitatm             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitatp             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitatu             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ habitatw             : num  0 0 0 0 0 0 0 0 0 0 ...
 $ seasons              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ seasonu              : num  0 1 0 0 0 1 0 1 0 0 ...
 $ seasonw              : num  1 0 1 1 1 0 1 0 0 1 ...

library(caTools)
set.seed(18)

split <- sample.split(mushroom_num$class, SplitRatio = 0.8)
training_set_num <- subset(mushroom_num, split == TRUE)
test_set_num <- subset(mushroom_num, split == FALSE)

table(training_set_num$class)

    0     1 
21745 27110 

table(test_set_num$class)

   0    1 
5436 6778 

Para comenzar a aplicar el clasificador del vecino más cercano, haremos uso de la función knn() tomando como valor de k=5.

library(class)
set.seed(18)
pred_knn <- knn(train = training_set_num[, -1], test = test_set_num[, -1], cl = training_set_num$class, k = 5)

Una vez tenemos los resultados de las predicciones, cosntruimos la matriz de confusión.

summary(pred_knn)

   0    1 
5429 6785 

confusionM <- table(test_set_num$class, pred_knn)
confusionM

   pred_knn
       0    1
  0 5429    7
  1    0 6778

accuracy_knn <- sum(diag(confusionM)) / sum(confusionM)
accuracy_knn

[1] 0.9994269

Como resultado de la matriz de confusión, tenemos una precisión del 99%, resultando ser mejor que en la regresión logística.

library(ROCR)
pred_knn_roc <- prediction(as.numeric(pred_knn), as.numeric(test_set_num$class))
perf_knn_roc <- performance(pred_knn_roc, "tpr", "fpr")
perf_knn_auc <- performance(pred_knn_roc, "auc")

print(perf_knn_auc@y.values[[1]])

[1] 0.9993561

plot(perf_knn_roc)

Una vez obtenemos la curva ROC y el resultado del área de debajo de la misma, observamos que el resultado es muy positivo, ya que gráficamente se aproxima a la vertical.

Clasificación con Árbol de Decisión

Comenzamos aplicando el árbol de decisión mediante la función rpart() como se muestra a continuación.

library(rpart)

Attaching package: ‘rpart’

The following object is masked from ‘package:dendextend’:

    prune

set.seed(18)
dt_classiffier <- rpart(class ~ ., data = training_set)

library(rattle)

Loading required package: bitops
Rattle: A free graphical interface for data science with R.
Version 5.5.1 Copyright (c) 2006-2021 Togaware Pty Ltd.
Type 'rattle()' to shake, rattle, and roll your data.

tree2 <- rpart(class ~ ., training_set,
  method = "class",
  control = rpart.control(cp = 0.00001)
)
fancyRpartPlot(tree2)

pruned <- prune(tree2, cp = 0.01)
fancyRpartPlot(pruned)

Construimos la matriz de confusión al igual que en los casos anteriores.

pred_dt <- predict(dt_classiffier, newdata = test_set, type = "class")

confusionM <- table(test_set$class, pred_dt)
confusionM

   pred_dt
       e    p
  e 4257 1179
  p  804 5974

accuracy_dt <- sum(diag(confusionM)) / sum(confusionM)
accuracy_dt

[1] 0.8376453

Una vez tenemos el resultado para el valor de precisión en la predicción una vez aplicado el árbol de decisión, 83%, pasamos a aplicar de forma visual la construcción de la gráfica de la curva ROC y el cálculo del AUC. En este caso, continúa siendo mejor resultado que el obtenido mediante la regresión logística, pero quedando por debajo que el resultado obtenido al aplicar KNN.

library(ROCR)
pred_dt_roc <- prediction(as.numeric(pred_dt), as.numeric(test_set$class))
perf_dt_roc <- performance(pred_dt_roc, "tpr", "fpr")
perf_dt_auc <- performance(pred_dt_roc, "auc")

print(perf_dt_auc@y.values[[1]])

[1] 0.8322468

plot(perf_dt_roc)

Random Forests classiffier

Comenzamos aplicando el clasificador de Random Forest mediante la función llamada de la misma forma; es decir, randomForest(). En esta función, el valor correspondiente al parámetro ntree indica la cantidad de árboles de decisión que formarán parte del clasificador.

library(randomForest)

randomForest 4.7-1.1
Type rfNews() to see new features/changes/bug fixes.

Attaching package: ‘randomForest’

The following object is masked from ‘package:rattle’:

    importance

The following object is masked from ‘package:dplyr’:

    combine

The following object is masked from ‘package:ggplot2’:

    margin

set.seed(18)
rf_classiffier <- randomForest(class ~ ., data = training_set, ntree = 250)

Realizamos las predicciones sobre el conjunto de datos y construimos la matriz de confusión.

pred_rf <- predict(rf_classiffier, newdata = test_set, type = "class")

confusionM <- table(test_set$class, pred_rf)
confusionM

   pred_rf
       e    p
  e 5431    5
  p    0 6778

accuracy_rf <- sum(diag(confusionM)) / sum(confusionM)
accuracy_rf

[1] 0.9995906

De forma similar al clasificador KNN, obtenemos un valor de precisión del 99%. Definimos la curva ROC para este caso y realizamos el cálculo del área bajo la curva, los cuales nos simbolizan su altísima precisión.

library(ROCR)
pred_rf_roc <- prediction(as.numeric(pred_rf), as.numeric(test_set$class))
perf_rf_roc <- performance(pred_rf_roc, "tpr", "fpr")
perf_rf_auc <- performance(pred_rf_roc, "auc")

print(perf_rf_auc@y.values[[1]])

[1] 0.9995401

plot(perf_rf_roc)

Kernel SVM Classifier

Comenzamos a aplicar el clasificador Máquina de Soporte Vectorial haciendo uso de la función svm(). En esta función, los valores de los parámetros type y kernel hacen referencia al tipo de clasificador; es decir, que el kernel es de tipo radial y gaussiano.

library(e1071)
set.seed(18)
svm_classiffier <- svm(class ~ .,
  data = training_set,
  type = "C-classification", kernel = "radial"
)

La predicción y la matríz de confusión son entonces las mostradas a continuación.

pred_svm <- predict(svm_classiffier, newdata = test_set, type = "class")

confusionM <- table(test_set$class, pred_svm)
confusionM

   pred_svm
       e    p
  e 5129  307
  p  240 6538

accuracy_svm <- sum(diag(confusionM)) / sum(confusionM)
accuracy_svm

[1] 0.9552153

Como se puede observar el valor de la predicción en este caso resulta ser del 95%. Para finalizar, construimos la curva ROC correspondiente en este caso

library(ROCR)
pred_svm_roc <- prediction(as.numeric(pred_svm), as.numeric(test_set$class))
perf_svm_roc <- performance(pred_svm_roc, "tpr", "fpr")
perf_svm_auc <- performance(pred_svm_roc, "auc")

print(perf_svm_auc@y.values[[1]])

[1] 0.954058

plot(perf_svm_roc)

Conclusiones

Una vez aplicados los cinco distintos métodos de clasificación sobre nuestro dataset tras realizar el preprocesamiento, podemos concluir diciendo que los clasificadores de KNN y Random Forest son los que poseen mayor precisión y, por tanto, menor error de predicción. Sin embargo, con los resultados obtenidos por parte de ambos podemos detectar un posible sobreajuste, no beneficiando al modelo. Por esta razón, creemos que resulta más beneficioso sacrificar parte del valor de precisión, teniendo en cuenta algunos valores de falsos positivos y falsos negativos, como ocurre en el caso del clasificador de la Máquina de Soporte Vectorial, aprovechando así su capacidad de mayor generalización.

accuracy_comp <- matrix(c(accuracy_rl, accuracy_knn, accuracy_dt, accuracy_rf, accuracy_svm), ncol = 5)
barplot(accuracy_comp)

barplot(accuracy_comp,
  main = "Accuracy Comparison",
  xlab = "Accuracy (%)",
  ylab = "Method",
  names.arg = c("RL", "K-NN", "DT", "RF", "SVM"),
  col = "darkred"
)

perf_auc <- matrix(c(perf_rl_auc@y.values[[1]], perf_knn_auc@y.values[[1]], perf_dt_auc@y.values[[1]], perf_rf_auc@y.values[[1]], perf_svm_auc@y.values[[1]]), ncol = 5)
barplot(perf_auc)

barplot(perf_auc,
  main = "AUC Comparison",
  xlab = "AUC (%)",
  ylab = "Method",
  names.arg = c("RL", "K-NN", "DT", "RF", "SVM"),
  col = "darkred"
)

Análisis aprendizaje no supervisado

En este apartado se analizará el dataset a través de algoritmos de aprendizaje no supervisado. En concreto, se probarán los algoritmos k-means y clustering jerárquico. Para ambos algoritmos, se seguirá el siguiente esquema:

1. Representar la distribución inicial de los datos.
2. Determinar el número de clusters óptimo.
3. Representar la distribución de los datos en función del número de clusters.
4. Calcular el promedio de las variables en función del cluster al que pertenecen.
5. Calcular el accuracy del algoritmo.

Para poder trabajar con algoritmos no supervisados será necesario que las variables sean numéricas. Para ello, se eliminarán las variables categóricas del dataset.

numerical_columns <- mushroom[, numerical_features]

K-means

En primer lugar, representaremos de forma gráfica la distribución inicial de los datos.

mushroom <- dummyVars(" ~ .", data = mushroom, fullRank = TRUE) %>% predict(mushroom)
mushroom <- as.data.frame(mushroom)

df <- as.data.frame(numerical_columns)

plot_ly(df,
  x = ~cap.diameter, y = ~stem.height,
  z = ~stem.width
) %>%
  add_markers(size = 1.5)

Para aplicar el algoritmo kmeans se utilizará la función kmeans() de la librería cluster. Será necesario determinar el número de clusters óptimo. Para ello, se utilizará la función “kmeans” con diferentes valores de “centers” y se calculará la suma de cuadrados internos (within groups sum of squares) para cada valor de “centers”. A continuación, se representará la suma de cuadrados internos vs. número de clusters.

wss_per_k <- 0
for (i in 1:10) {
  kmeans_aux <- kmeans(numerical_columns, center = i, nstar = 20)
  wss_per_k[i] <- kmeans_aux$tot.withinss
}

Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)

par(mfrow = c(1, 1))
plot(1:10, wss_per_k,
  type = "b",
  xlab = "Number of clusters",
  ylab = "WSS",
)

Como se puede observar en la gráfica anterior, la suma de cuadrados internos disminuye a medida que aumenta el número de clusters. Sin embargo, a partir de 2 clusters, la disminución de la suma de cuadrados internos es muy pequeña. Por lo tanto, se decide utilizar 2 clusters. Comprobamos que tiene sentido utilizar 2 clusters, ya que conocemos que el dataset es binario.

Generamos el modelo de k-means con 2 clusters.

km_model <- kmeans(df, center = 2, nstar = 20)

Representamos la distribución de los datos en función de los clusters obtenidos.

df$cluster <- factor(km_model$cluster)

plot_ly(df,
  x = ~cap.diameter, y = ~stem.height,
  z = ~stem.width, color = ~cluster
) %>%
  add_markers(size = 1.5)

Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels

Se puede observar que los champiñones de menor tamaño (incluyendo diámetro, altura y anchura) pertenecen al cluster 2 y los de mayor tamaño pertenecen al cluster 1.

A continuación, calcularemos el valor promedio de las variables para cada cluster generado con el modelo de k-means.

grouped_mushroom <- df %>%
  group_by(cluster) %>%
  summarise(
    mean_cap_diameter = mean(cap.diameter),
    mean_stem_height = mean(stem.height),
    mean_stem_width = mean(stem.width)
  )

grouped_mushroom

A partir de este momento, hemos decidido modificar el dataset actual debido a que para aplicar técnicas como “silhouette” o “dendrogram” (para el caso de clustering jerárquico) es necesario que el dataset sea de menor tamaño. Para mantener la proporción de los datos de cada clase, se reducirá el dataset haciendo uso de createDataPartition, manteniendo únicamente un 1% de los datos iniciales para el análisis.

set.seed(42)
split <- createDataPartition(mushroom$class, p = 0.01)
smaller_df <- mushroom[split$Resample1, ]

Comprobamos que la proporción de datos de cada clase se mantiene al hacer la partición.

initial_class_prop <- table(mushroom$class) / nrow(mushroom)
smaller_class_prop <- table(smaller_df$class) / nrow(smaller_df)

print(initial_class_prop)

        0         1 
0.4450867 0.5549133 

print(smaller_class_prop)

        0         1 
0.4386252 0.5613748 

Ya podemos trabajar con el dataset reducido. Comprobamos que el dataset ahora cuenta con 611 ejemplos y 3 atributos.

smaller_df <- smaller_df[, numerical_features]
dim(smaller_df)

[1] 611   3

En primer lugar, visualizaremos la distribución inicial de los datos a través de una gráfica 3D.

plot_ly(smaller_df,
  x = ~cap.diameter, y = ~stem.height,
  z = ~stem.width
) %>%
  add_markers(size = 1.5)

A continuación, vamos a estudiar cuál sería el número óptimo de clusters para el dataset reducido haciendo uso de la medida de bondad interna “silhouette”. Para ello, utilizaremos la función fviz_nbclust de factoextra. Silhouette es una medida que sirve para validar el número de clusters. Se calcula como la diferencia entre la distancia media de un punto a los puntos de su propio cluster y la distancia media de un punto a los puntos de su cluster más cercano.

fviz_nbclust(smaller_df, FUNcluster = kmeans, method = "silhouette")

Según la gráfica, podemos afirmar que el número óptimo de clusters es 2.

Por último, visualizamos la distribución de los datos en función de los clusters obtenidos.

km_sm_model <- kmeans(smaller_df, center = 2, nstart = 20)
cluster <- factor(km_sm_model$cluster)

plot_ly(smaller_df,
  x = ~cap.diameter, y = ~stem.height,
  z = ~stem.width, color = ~cluster
) %>%
  add_markers(size = 1.5)

Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels

Por último, calculamos el valor promedio de las variables para cada cluster generado con el modelo de k-means.

grouped_sm_mushroom <- smaller_df %>%
  mutate(cluster = cluster) %>%
  group_by(cluster) %>%
  summarise(
    mean_cap_diameter = mean(cap.diameter),
    mean_stem_height = mean(stem.height),
    mean_stem_width = mean(stem.width)
  )
grouped_sm_mushroom

Clustering jerárquico

A continuación, vamos a aplicar el algoritmo de clustering jerárquico a nuestro dataset reducido. Para ello, utilizaremos la función hclust. Primero, calculamos la distancia entre los puntos del dataset.

distance <- dist(smaller_df)
hc_model <- hclust(distance)

Representamos el dendrograma para visualizar la distribución de los datos en función de los clusters obtenidos.

dend_modelo <- as.dendrogram(hc_model)
plot(dend_modelo)

Hasta ahora, hemos obtenido la jerarquía de los datos, pero lo que realmente nos interesa es la clasificación de los datos en función de los clusters. Cortaremos el dendrograma en un punto que nos interese para obtener los clusters. En este caso, hemos decidido cortar el dendrograma en 90 para obtener una visualización del dendograma cortado.

cut <- 0.9

dend_modelo %>%
  color_branches(h = cut) %>%
  color_labels(h = cut) %>%
  plot()

Para obtener el número óptimo de cluster, haremos uso de la medida interna de bondad silhouette. Para ello, utilizaremos la función fviz_nbclust de factoextra.

fviz_nbclust(smaller_df, FUNcluster = hcut, method = "silhouette")

Comprobamos que en este caso, el número óptimo de clusters podría ser 2 o 3, ya que el valor de silhouette es muy similar para ambos casos. En este caso, hemos decidido utilizar 2 clusters para poder comparar posteriormente los resultados con los obtenidos con el algoritmo de k-means.

Calculamos la agrupación del modelo en función del número de clusters que hemos decidido utilizar. Además, calculamos el promedio de los datos de cada cluster para ver si podemos sacar alguna conclusión.

jq_cluster <- cutree(hc_model, k = 2)

grouped_mushroom <- smaller_df %>%
  mutate(cluster = jq_cluster) %>%
  group_by(cluster) %>%
  summarise_all(mean)
grouped_mushroom

Visualizamos la agrupación de los datos en función de los clusters obtenidos.

jq_cluster <- factor(jq_cluster)

plot_ly(smaller_df,
  x = ~cap.diameter, y = ~stem.height,
  z = ~stem.width,
  color = ~jq_cluster
) %>%
  add_markers(size = 1.5)

Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels

Con el objetivo de comparar los resultados obtenidos en los dos algoritmos, vamos a calcular el rendimiento de cada uno de ellos, haciendo uso del accuracy como medida de bondad externa.

En primer lugar, calculamos el accuracy del modelo de k-means. Supondremos que la clase 1 es la clase “e” y la clase 2 es la clase “p”. Para ello, obtenemos las clases reales y las clases predichas, y calculamos el accuracy.

Volvemos a obtener el dataset reducido para poder tener las clases reales.

smaller_df <- mushroom[split$Resample1, ]

real_classes <- ifelse(smaller_df$class == "e", 1, 2)
predicted_classes <- km_sm_model$cluster
predicted_classes <- as.numeric(predicted_classes)

accuracy <- sum(real_classes == predicted_classes) / length(real_classes)
print(accuracy)

[1] 0.700491

Hacemos lo mismo con el modelo de clustering jerárquico, pero en este caso, supondremos que la clase 1 es la clase “p” y la clase 2 es la clase “e”.

real_classes <- ifelse(smaller_df$class == "e", 2, 1)
predicted_classes <- as.numeric(jq_cluster)

accuracy <- sum(real_classes == predicted_classes) / length(real_classes)
print(accuracy)

[1] 0.9885434

Tras comparar los resultados obtenidos en los dos algoritmos, podemos afirmar que el modelo de clustering jerárquico ha obtenido un accuracy mayor para este dataset.